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ABSTRACT 

Sub milli-arcsecond imaging in the visible band will provide a new perspective in stellar astrophysics. Even 
though stellar intensity interferometry was abandoned more than 40 years ago, it is capable of imaging and 
thus accomplishing more than the measurement of stellar diameters as was previously thought. Various phase 
retrieval techniques can be used to reconstruct actual images provided a sufficient coverage of the interferometric 
plane is available. Planned large arrays of Air Cherenkov telescopes will provide thousands of simultaneously 
available baselines ranging from a few tens of meters to over a kilometer, thus making imaging possible with 
unprecedented angular resolution. Here we investigate the imaging capabilities of arrays such as CTA or AGIS 
used as Stellar Intensity Interferometry receivers. The study makes use of simulated data as could realistically 
be obtained from these arrays. A Cauchy-Riemann based phase recovery allows the reconstruction of images 
which can be compared to the pristine image for which the data were simulated. This is first done for uniform 
disk stars with different radii and corresponding to various exposure times, and we find that the uncertainty 
in reconstructing radii is a few percent after a few hours of exposure time. Finally, more complex images are 
considered, showing that imaging at the sub-milli-arc-second scale is possible. 

1. INTRODUCTION 

There has been a recent interest in the revival of Stellar Intensity Interferometry (SII) due to the excellent 
baseline coverage of planned air Cherenkov telescope arraysP^ This interest has lead to developments in 
instrumentation, experimentation, and simulations of the capabilities of this technique. Various analog and 
digital correlator technologies^ are being implemented by LeBohec et. al.J^and cross correlation of streams of 
photons with nanosecond scale resolution has already been achieved. The suitability of various proposed array 
configurations is being evaluated by Jensen et alP to understand their different sensitivities for interferometric 
imaging before final choices of the array layouts are made. Image reconstruction algorithms such as the one 
suggested by Holmes et. alP have opened the possibility of imaging and have become interesting subjects 
in their own right. In this paper we will focus on the imaging capabilities and limitations of air Cherenkov 
telescope arrays used as high angular resolution intensity interferometers. 

High angular resolution astronomy in the optical range will open a whole new field of exploration. The 
possibility of viewing stars as extended objects will enable the testing of many current astrophysical models, 
and the knowledge acquired will have consequences on many fields related to stellar astrophysics. As a first 
example consider the measurement of stellar diameters, which can be performed to an accuracy of a few percent 
with the methods discussed in this paper (see section |4~T]) . When measuring diameters at different wavelengths, 
we can learn about the behaviour of the optical depth as a function of the radius of a particular star. 7 This 
type of measurement becomes particularly useful at shorter wavelengths (A ~ 400nm) than those feasible with 
conventional amplitude (Michelson) interferometry. Another interesting science case is the study of fast rotating 
stars, for which we can measure oblateness, pole brightening, and disk formation. There is also the study of 
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interacting binaries. An actual image of an interacting binary will not only aid in the determination of the or- 
bital parameters to a high degree of accuracy, but also in the study and imaging of mass transfer and accretion. 
Although not rigorously discussed in this paper, imaging mass transfer will improve our understanding of late 
stellar evolution, i.e. the formation of type II supernovae, which are of great importance in cosmologyp' and our 
understanding of the formation of compact objects. The number of interesting science cases and astrophysical 
targets is overwhelmingly large, and a detailed discussion is given by Dravins et alP 



We propose the use of Intensity Interferometry for high angular resolution astronomy (see Holder et. alP 
for more details). This technique was introduced in the 1950's by R. Hanbury BrowrP^Ell and implemented in 
the 1960's with the Narrabri Stellar Intensity Interferometer pi accomplishing the measurement of over thirty 
stellar diameters. The use of planned air Cherenkov telescope arrays poses a unique opportunity to revive 
Intensity Interferometry. With hundreds of telescopes separated by up to ~ lkm, it will be possible to have an 
unprecedented coverage of the Fourier plane and thus achieving sub-milliarcsecond resolution. 



The most important difference between SII and amplitude interferometry is that SII relies on the correlation 
between the low frequency intensity fluctuations and so does not rely on the relative phase of the individual 
waves received at different telescopes. Intensity interferometry measures the squared modulus of the complex 
degree of coherence I7I 2 . 

M a = <A/ ' Af '-> (i) 

Here < > is the time average of the intensity received at a particular telescope i, and Ali refers to the 
low frequency intensity fluctuations received at telescope i. Intensity interferometry has several advantages and 
disadvantages when compared to amplitude interferometry. The main advantages are that it is insensitive to 
atmospheric turbulence and that it does not require high optical precisioifj]. The complex degree of coherence 
is proportional to the Fourier transform of the object in the sky (Van Cittert-Zernike theorem), and since we 
measure the modulus squared of 7, the main disadvantage is that the phase of the Fourier transform is lost 
in the measurement process, posing difficulties in recovering actual images from magnitude information only. 
In addition to imaging difficulties, measuring a second order effect also results in sensitivity issues, 1 which can 
be dealt with by using large light collection areas and exposure times. As for the imaging limitations, several 
phase retrieval techniques exist, and we will implement a two dimensional version of the one dimensional ap- 
proach introduced by Holmes & Belen'kii. 6 It is important to note that our results pertain to a single phase 
recovery algorithm,^ and a comparison to other algorithms is currently being investigated.^ Once a sufficient 
coverage of the Fourier plane is available, and phase recovery is performed, a study on imaging capabilities 
can be performed. Here we will first concentrate on the study of the uncertainty when reconstructing disk-like 
stars. Then a less exhaustive analysis on the capabilities is performed for more complicated images such as 
oblate rotators, binary stars and stars with bright & dark features. 



The outline of the paper is the following: First we will briefly discuss the phase retrieval technique. Then 
we will discuss the simulation of our data and how it will be fitted to an analytic function so that the phase 
retrieval method can be applied. Finally we will discuss the capabilities for imaging disk-like stars, binary stars 
and more complicated objects. 



2. PHASE RECONSTRUCTION 

The objective of phase reconstruction is to recover the phase of the Fourier transform of the image from magni- 
tude information onlyP The resulting image is then reconstructed up to an arbitrary translation and reflection. 
It is simpler to first understand phase retrieval in one dimension and then generalize to two dimensions . One 
possible route towards phase retrieval starts by first approximating the continuous Fourier transform I(x) by 

T For a more detailed discussion on the advantages of Intensity Interferometry see Hanbury Brown 12 



a discrete one (I(mAx) — Q(j A9)e^ mk ° AxAe , where 0(6) is the image in the sky and ko is the usual wave 
vector). Then the discrete Fourier transform can be expressed as a magnitude times a phasor (I(z) = R(z)e 1 ^^ 
where z = e " nk o^A9 j g com pi ex ). The most important step is then to apply the theory of analytic functions 
i.e. the Cauchy-Riemann equation!!. These relate the phase $ and the log- magnitude lni? along the real or 
imaginary axes. One can show by using the Cauchy-Riemann equations, that the phase differences along the 
radial direction in the complex plan^l are directly related to the differences in the logarithm of the magnitude 
(see Holmes & Belen'kiP for more details) , so that integrating the Cauchy-Riemann equations directly does not 
immediately solve for the phase. In other words, phase differences along the purely real or imaginary axes are 
not available directly from the data. 

Since z, the independent variable of the Fourier transform (z = e lmk o^xA9^ ^ ag moc [ ums e q ua l to 1, the 
phase differences that we seek lie along the unit circle in the complex plane. Consequently, the procedure to 
find the phase consists in first assuming a plausible solution form, then taking differences in the radial direction 
of the complex plane, and finally fitting the data to the radial differences of the assumed solution. A general 
form of the phase can be postulated by noting that the phase is a solution of the Laplace equation in the 
complex plane (applying the Laplacian operator on the phase and using the Cauchy-Riemman equations yields 
zero) . Since the phase differences are known along the radial direction in the complex plane we can take radial 
differences of the general solution and then fit the log-magnitude differences (available from the data) to the 
radial differences of the general solution. 

One can think of this one-dimensional reconstruction as a the phase estimation along a single slice in the 
Fourier plane. A generalization to two dimensions can be made by doing the same procedure for several slices. 
The direction of the slices is arbitrary, however for simplicity we reconstruct the phase along horizontal or 
vertical slices in the Fourier plane, and noting that one can relate all slices with a single orthogonal slice, i.e. 
once the phase at the origin is set to zero, the single orthogonal slice sets the initial values for the rest of the 
slices. The resulting reconstructed phase will be arbitrary up to a constant and a linear term, which corresponds 
to a translation. It should be noted that the above solution approach gave reasonably good results. However, 
it is not the only possible approach. We have also investigated Gerchberg-Saxton phase retrieval, Generalized 
Expectation Maximization, and other variants of the Cauchy-Riemann approach.^ 

3. PROCEDURE 

Having briefly discussed the phase reconstruction algorithm, our basic procedure for recovering images is the 
following: First we simulate realistic data as would be obtained from an air Cherenkov telescope array such 
as CTA or AGIS. Once simulated data are available, they are fitted to an analytic function so that the phase 
recovery algorithm can be applied in a more straightforward way. Finally, once the phase is recovered, the 
inverse Fourier transform will provide us with a reconstructed image. Some details concerning the simulations 
and data fitting approach now follow. 

3.1 Simulation of realistic data and array design used 

The simulation of the the data that may be produced by an array of telescopes starts from a pristine image, 
generally 2048 x 2048 pixels with an arbitrary dynamic range. The squared magnitude of the Fourier transform 
of this image is obtained by means of an FFT algorithm and it is normalized so its maximum at zero baseline 
is equal to one. With a wavelength A = 400 nm, and the full scale of the image typically set to 10 mas, this 
provides a value for the expected degree of coherence on a square grid with a pitch of ~ 8.2 m extending over 
a ~ 16.8 km x 16.8 km area. 

The squared Fourier magnitude map is sampled by the set of pairs available in the simulated array. Sim- 
ulations presented here have been obtained with an array of N — 97 telescopes, each with a light collecting 

'''The C-R equations can be applied because I(z) is a polynomial in z, where z = e 1 ^ . 

^If £ is the real axis and ip is the imaginary axis, then a difference along the radial direction is A£ + iAi/>. 
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Figure 1. Array configuration used for our analysis. 



area of 100m 2 and a quantum efficiency a — 0.30 resulting in an effective area of 30 m 2 . The telescopes are 
distributed in the field according to an early design of the CTA array shown in figure [2] Such an array provides 
a coverage of the interferometric (u,v) plane with N(N — l)/2 = 4656 baselines many of which are redundant. 
The baselines are shown in figure [5] The degree of coherence recorded by each baseline is obtained from a linear 
interpolation between the closest four points in the Fourier magnitude table. 

The data recorded by a real array would be affected by the diurnal motion of the observed star which affects 
the effective baselines by projection. The average correlation must then be recorded for each baseline at time 
intervals short enough for the baseline change to be negligible. In the first simulation study reported here we 
have decided to avoid this complication and simulate data that would result from the observation of fixed stars 
at zenith so the effective baselines used are those shown in figure [5] without any further projective distortion. 
The implications of this simplification choice are a less uniform sampling of the (u,v) plane compensated by 
smaller error bars on the degree of coherence from each baseline record. These two effects essentially cancel each 
other as long as small scale features in the (u,v) plane are not central to the analysis. The benefits from the 
simplification is a reduction in the volume of data to handle (each simulation produces a single record for each 
pair of telescopes) and eliminates further arbitrary parameters (such as the site latitude, celestial declination, 
range of hour angles and time interval between recordings). 

Once the degree of correlation within each baseline has been obtained, Gaussian noise is added. The 
Gaussian nature of the noise was tested with detailed simulation of a pair of photo-multiplier tube signals 
corresponding to a random stream of photons. The time integrated product of the two traces was Gauss 
distributed. The magnitude my of the star is used to compute a spectral density (m~ 2 s~ 1 Hz~ 1 ) according 
to n = 5 x I0~5 x 2.5~" lv \ The standard deviation of the Gaussian noise added to the pair of telescopes 
(i, j) is calculated as a = n^J Ai ■ Aj ■ Af • At/2 where At is the effective light collection area of the i th tele- 
scope, Af is the signal band-width and At is the observation duration. For simulations presented in this paper 
A/ = 200MHz which is a realistic choice when considerations on air Cherenkov telescope optics and electronics 
are taken into account (See Holder & LeBohecPand LeBohec et alPJ). 
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Figure 2. Available baselines for the array design used in this study. 
3.2 Fitting the data to an analytic function 

The phase reconstruction algorithm is greatly simplified when data are known on a square grid rather than in 
a 'randomly' sampled way as is directly available from observations. This is because sampling data on a fine 
square grid enables an easier estimation of the derivatives of the log-magnitude. 



Assuming that data f{xi) are known at positions Xi 
function that minimizes the following \ 2 
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Here, the ' s arc the coefficients of the basis functions that we want to use to fit our simulated data. 
Any complete basis will suffice in theory, however it is more appropriate to choose a set of basis functions 
that tend to zero at infinity. The reason for this requirement of our basis functions is so that data are more 
realistically fitted in regions were there is not much data available. For the case of the CTA array design that 
we used, we noticed that there is less data at baselines greater than 600 m (see figure [5]) . We found that basis 
functions that meet this requirement are the solutions to the two dimensional quantum-mechanical harmonic 
oscillator, i.e. Hermite polynomials with Gaussian envelopes. These also turn out to be convenient because 
they are eigenfunctions of the Fourier transform operator. 

Now we can turn this problem into a linear system by taking derivatives with respect to the unknowns a^. 
Our data fitting typically starts with a small number of basis elements, then wc check to sec if a certain reduced 
X 2 is met (in our case, we chose an acceptable reduced \ 2 to be 1.5), if this condition is not met, then the 
number of basis elements is increased iteratively until it does. 



4. IMAGING CAPABILITIES 

We may start by quantifying the resolution of an air Cherenkov telescope array such as the one shown in figure 
[2] by recalling that quantities that are related to each other by a Fourier transform obey an uncertainty relation. 
For an order of magnitude estimation it suffices to relate the size of the array to the maximum resolution by 



(3) 



With a kilometer size array (Ax ~ lkm), and a wavelength of A ~ 400nm we obtain a resolution of 
AO ~ 0.1 mas. On the other hand, the largest objects observable with an array whose inter-telescope sepa- 
ration is of the order of Ax ~ 50m is AO ~ 1 mas. These order of magnitude considerations will be taken 
into account when performing simulations and image reconstructions, i.e. the minimum and maximum size of 
pristine images will not go far beyond these limits. 

We tested the imaging capabilities for simple objects, namely uniform disk-like stars, oblate rotating stars, 
binaries, and more complex images. First we will concentrate on the capabilities and limitations for recon- 
structing uniform disk-like stars. We will show that such a preliminary analysis reveals more precisely, when 
compared to the previous estimate, the sizes of objects that can be observed. Even though using a Cauchy- 
Riemann based approach to recover images might not be the most efficient way to measure stellar radii, such 
a study will start to quantify the abilities of measuring other scale parameters in more complicated images 
(oblateness, distance between binary components, etc.). 

4.1 Uniform disks 

Simulated data were generated corresponding to uniform disk stars of a particular brightness. We set the 
brightness to magnitude 6 after noting that an error of a few percent in the simulated data can be achieved in 
a few hours. Also, 6 th magnitude stars are appropriate since they roughly correspond to the upper limit for 
most of the interesting astrophysical targets found by Dravins et alP 

It is interesting to first study the uncertainty for a particular exposure time and a brightness corresponding 
to 6 th magnitude. We simulated data corresponding to uniform disks of random radii up to 1 mas for 50 hours 
of exposure time. An example of such a reconstruction is shown in figure [3Jd, where the brightness is shown 
in arbitrary units between and 1. A first look at the example reconstruction reveals that the edge of the 
reconstructed disk is not sharp, so a threshold in the brightness was applied for for measuring the radius. The 
radius was measured by counting pixels above a threshold and noting that the area of the disk is proportional 
to the number of pixels passing the threshold. After experimenting for different radii, we chose the threshold for 
measuring the radius to be 0.2. We can now compare the simulated and reconstructed radii as shown in figure 
|3t, where each point in the figure corresponds to an individual simulation (including noise) and reconstruction. 

Figure |3Ji clearly shows that stellar diameters ranging from 0.05 mas to 0.5 mas can be measured with 
uncertainties smaller than 5%. The 'tail' seen in the bottom left of figure |3Ji shows the smallest measurable 
radius, so we take this radius (~ 0.03 mas) to be the point spread function (PSF) of the array. The uncertainty 
shown in the sub-figure in figure [3] was estimated after running noisy simulations many times as is shown in 
figure U] 

It can be seen from figure [3^, that the uncertainty increases roughly linearly as a function of the pristine 
(simulated) radius. This can be understood from the following argument: As the pristine radius decreases, the 
distance to the first zero in the correlation increases as ~ r, so the number of telescopes contained within the 
Airy disk increases as ~ r 2 . Consequently, decreasing the pristine radius is equivalent to increasing the number 
of independent measurements by a factor of ~ r 2 at most (at least by a factor of r). Since the uncertainty will 
decrease as the square root of the number of independent measurements, the error will decrease as ~ r at most. 
For radii above 0.6 mas, there are simply not enough baselines to constrain the Fourier plane information for 
image reconstruction. For radii greater than 0.6 mas, the size of the Airy disk is of the order of 100 m, and this 
results in having less than 100 independent measurements as can be seen in figure [2] Also shown in figure 0] is 
the percent error as a function of time for two different radii, where it can be seen that a percent error of less 
than 5% is achieved after only a few hours. 

5. COMPLEX IMAGES 

Our algorithm has also been tested on more complex images such as oblate rotating stars, binary stars, and 
stars with brighter or darker regions. Since we have not developed a proper tool for quantitatively comparing 
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Figure 3. a) Simulated vs. Reconstructed radii for magnitude 6 stars with 50 hours of observation time. The top sub- 
figure shows the uncertainty for a 0.2 mas measurement. The bottom sub-figure shows the residual (Reconstructed-Real) 
along with the uncertainty in the radius, b) Example of a reconstructed uniform disk of radius 0.2 mas. Also shown is 
a slice of the reconstructed image (solid line) compared to a slice of the pristine image (dashed line). 




Figure 4. a) Histogram of real radius minus reconstructed radius for 50 hours of exposure time on a 6 th magnitude star 
of 0.2 mas radius, b) Percent error as a function of time for several reconstructed radii. 
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Figure 5. a) Simulated and reconstructed oblate rotator of magnitude 3 and 10 hours of observation time, b) Simulated 
and reconstructed binary of magnitude 6 and 12 hours of observation time. 



simulated and reconstructed imagefl we will only show a few representative examples of what the algorithm is 
capable of. Figure [5£i shows an example reconstruction of an oblate rotating star of magnitude 3 and 10 hours 
of observation time, the semi-major axis and semi-minor axis of the pristine ellipse are 0.2 mas and 0.12 mas 
respectively. Also shown in this figure, and all subsequent ones, is the pristine image convolved with the point 
spread function of the array (see section l4~Tj) . Reconstructed noise of less than 10% can be seen in the figure. 
These noise fluctuations are not so much a consequence of noise in the simulated data, but of the reconstruction 
algorithm itself, and the preferential direction of this reconstructed noise (along the vertical direction) is due 
to our choice of the slice direction for the phase reconstruction, i.e. the phase was reconstructed by taking 
horizontal one-dimensional slices of the magnitude, and then related to each other with a single vertical slice. 
As for the structure (bumps) within the star in figure Ek, these start to appear when either the star becomes 
bright or enough exposure time is supplied so that information other than the first lobe in the Fourier magni- 
tude is significant. In other words, when high frequency portions are visible in the Fourier magnitude, fictitious 
structure starts to become visible. This is most likely due to the fact that most of the high frequency infor- 
mation in the Fourier plane is used to reconstruct a dark background of several mas's with a central bright region. 

The case of a binary star is shown in figure Eb, corresponding to a magnitude 6 binary star for 12 hours 
of exposure time. The noise in the reconstructed image has the same origin as in the case of the oblate star. 
Although the inclination angle was well reconstructed in both cases in figureEl image reconstruction is degraded 
when the symmetry axis is neither the x or y axis. Again, this is due to the particular phase recovery method 
of taking horizontal or vertical slices and the degradation is significantly reduced when aligning one symmetry 
axis of the magnitude to our x or y axis. 

Having this symmetry consideration in mind, a slightly more complicated example is shown in figure EK, 
corresponding to a star obscured by a disk (of dust for example). A black streak in the pristine image repre- 
senting the obscuring disk is aligned with the x axis. The black streak can be easily seen in the reconstruction 
as well as the contour of the obscured star. This image becomes increasingly easier to reconstruct as the image 
becomes more and more symmetric, that is, as the black streak in the pristine image moves towards the center 



" Since the simulated and reconstructed objects can differ by translations and reflections, developing a tool than can 
accurately quantify the difference between simulated and reconstructed images is not trivial. 
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Figure 6. a) Simulated and reconstructed star obscured by a disk. This corresponds to magnitude 4 and 15 hours of 
observation time, b) Uniform disk of radius 0.2 mas with a dark spot of radius 0.05 mas. The simulation was done for 
magnitude 6 and 100 hours of observation time. 

of the staJl. As a final example, we considered the case of a dark spot in a star as shown in figure [6b- The 
simulation corresponds to a magnitude 6 star and 100 hours of observation. Even though the reconstructed 
image appears brighter in the bottom, and the location of the spot is more symmetric, the size of the spot 
(comparable to the PSF) is well reconstructed. 

6. CONCLUSION AND OUTLOOK 

We have shown that planned air Cherenkov telescope arrays have sufficient baselines to provide an excellent 
coverage of the (u,v) plane. By first simulating realistic data, we show that it is possible to achieve a signal to 
noise ratio of the order of 10 within a few hours for relatively faint stars (magnitude 6), and obviously higher S/N 
for brighter stars. Once data were simulated we also show that imaging is possible using a Cauchy-RiemanrP 
phase recovery technique. A study of the error propagation with disk-like stars reveals that the uncertainty in 
reconstructed radii is of a few percent. We also explored imaging capabilities for more complex images such 
as oblate rotating stars, binary stars, and stars with dark/bright regions yielding good results. A quantitative 
analysis of the reconstruction capabilities for complex images is still in progress. 

The array design could certainly be improved by including telescopes at shorter distances. This could sig- 
nificantly improve the size range of observable objects, in particular, we could observe objects of more than 1 
mas across at 400 nm. 



The analysis can also be improved by doing several things: If the pristine object has a symmetry axis, a 
first reconstruction can be made to find it, and then a second image reconstruction can be made to improve the 
first results. Something else worth implementing is a first reconstruction only constraining the low frequency 
components of the phase by using the low frequency part of the magnitude, then a second reconstruction could 
be performed, only dealing with the internal details of the image. 



"This can be understood by noting that the Fourier magnitude of a non symmetric binary becomes almost indistin- 
guishable from the one corresponding to a central bright star with two fainter companions at either side. 



To conclude this simulation phase, pristine images generated from astrophysical models should be gener- 
ated in order to identify how much of the astrophysical model can actually be constrained. This aspect of the 
simulation phase is currently under development. 
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